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Population extinction is a rare event which requires overcoming an effective barrier. We show that 
the extinction rate can be fragile: a small change in the system parameters leads to an exponentially 
strong change of the rate, with the barrier height depending on the parameters nonanalytically. 
General conditions of the fragility are established. The fragility is found in one of the best-known 
models of epidemiology, the SIS model. The analytical expressions are compared with simulations. 



PACS numbers: 87.23.Cc, 05.40.-a, 02.50.Ga 

Extinction of a population is of central interest for pop- 
ulation dynamics [J 0] • It results from a large fluctua- 
tion away from a steady state of the population. Such 
fluctuations are usually rare. They require an unlikely se- 
quence of elementary birth-death events or a large change 
in the fluctuating environment, or both, often visualized 
as overcoming a barrier. Much work has been done on 
extinction for various fluctuation mechanisms, and the 
extinction rates have been found for a number of models 
of population dynamics @, i, 0, @, 0, 0, S, M, d E3, d • 

There is a close similarity between population extinc- 
tion and a diverse group of physical phenomena which 
involve switching between coexisting states and range 
from switching in Josephson junctions and nanomagnets 
to chemical reactions and to protein folding. Both extinc- 
tion and switching are caused by large rare fluctuations. 
In many cases the rate of extinction (switching) W is ex- 
ponentially small, W oc exp(— Q) with Q » 1 [T3|. In 
particular, for systems close to thermal equilibrium the 
switching exponent Q is Q = R/ksT, where R is the free 
energy barrier and T is temperature fl5| . 

In this paper we show that, in spite of the aforemen- 
tioned similarity, population extinction displays a fea- 
ture that does not generally occur in switching between 
metastablc states. We find that the extinction rate is of- 
ten fragile. A small perturbation of the system can lead 
to an abrupt change of the rate exponent Q. If the per- 
turbation is proportional to a parameter /i, the value of 
Q for \i = is much larger than for /j — > 0. We find a gen- 
eral condition for the fragility to occur and illustrate the 
effect with a broadly used epidemiological model, the so 
called SIS model where there are present only susceptible 
(S) and infected (I) individuals Q. 

The difference between switching and extinction can 
be understood from Fig. [TJ For illustration purpose, sys- 
tems that can switch or display extinction are sketched as 
particles with dynamical variables x moving in a poten- 
tial U (x) , even though the actual system motion is gen- 
erally non-potential. In population dynamics, the com- 
ponents of x determine the size of different populations. 
In switching, if the system is initially near a stable state 
x^4, over the relaxation time t T there is formed a quasi- 
stationary current away from the basin of attraction to 
X.A- The current gives the switching rate [l5| . It goes over 
the saddle point [3, HE EE [13, EH and is divergence- free 



there, the probability distribution does not accumulate 
near the saddle point. 




FIG. 1: A sketch of switching between stable states (a) and 
extinction (b); x are the dynamical variables of the system. 
The stable states correspond to the minima of the effective 
potential U(x). In switching, a quasi-stationary probability 
current shown by arrows goes from the initially occupied to 
the initially empty state over the saddle point. In extinction, 
the probability current terminates at the extinction state and 
the occupation of this state linearly increases in time for t <C 
W , where W is the extinction rate. 

The extinction rate is also given by the quasi- 
stationary current away from the vicinity of the stable 
state, see Fig. [1] (b). The current goes to the hyper- 
plane xe — where population E goes extinct and, in a 
qualitative difference from the switching current, termi- 
nates there, because the population size may not become 
negative. Since the extinct population usually does not 
emerge again, fluctuations do not remove the system from 
the hypcrplane xe = 0. Therefore the probability distri- 
bution accumulates there, as seen in Fig.[5]Ja) below. For 
xe > the probability distribution is quasi-stationary 
on times t t < t < W' 1 . 

Because of the current discontinuity, near xe = 
the quasi-stationary distribution sharply varies with xe 
(exponentially sharply, see below). However, with re- 
spect to non-extinct populations, x^e, it generally has 
a smooth maximum, cf. Fig. [DJa). This is no longer 
true if the system has effective constraints, for example 
the total population is conserved. Population conserva- 
tion is a commonly used assumption in the SIS model 
and indeed the distribution over x^ E 
is sharp near xe = in the SIS model with constant 
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population. Once the constraint is lifted, for example, 
the population starts fluctuating, the shape of the dis- 
tribution changes as does also the extinction rate. The 
rate change occurs in an exponentially narrow parameter 
range and is exponentially large, which is a signature of 
the fragility. 

We consider extinction in a spatially uniform system of 
coupled populations (species). The system state is char- 
acterized by a vector X with integer components equal 
to the size of different populations. The quasi-continuous 
vector x = X/JV introduced above gives the populations 
size scaled by a large characteristic total population N. 
The dynamics is quite generally described by a master 
equation for the probability p(X), 

p(X) = [W(X - r; r)p(X - r) - W(X; r)p(X)] . (1) 

r 

Here, VF(X; r) is the rate of an elementary transition 
X — > X + r in which the populations change by r = 
(fx, r2, • • ■)• The condition that the system does not leave 
the extinction hyperplane has a form 

W(X; r) = for X E = 0, r E ± 0. (2) 

If fluctuations can be disregarded, from Eq. (JTJ) we ob- 
tain for average scaled populations x a mean-field equa- 
tion 

x = V rw(x;r), (3) 

L J r 

where u>(x;r) = W(X;r)/N is a characteristic transi- 
tion rate per individual. We assume that Eq. ([3]) has an 
asymptotically stable solution x^ and a stationary solu- 
tion xs that lies on the extinction hyperplane xe = 0, 
cf. Fig. DJb), and is asymptotically stable with respect 
to x^e but unstable for xe- For t r -c t <C W^ 1 the 
distribution p(X) peaks at X^ = Ntca- 

The exponent Q in the extinction rate W oc exp Q can 
be found by either solving the mean first passage time 
problem for reaching extinction J3j,jj, nfl or by calculat- 
ing the small-X B tail of p(X) @T 0, [lj, E3]. In both 
methods one looks for the optimal (most probable) fluc- 
tuation that leads to extinction, and the results coincide. 
Here we will study the quasi-stationary distribution. In 
a standard way, we seek the solution of Eq. ([T]) in the 
cikonal form, 

p(X) = exp[-Na(x)], a = -H(x,d x s), 

p) = V w(x; r) [exp(pr) - I] . (4) 

We took into account that, typically, |r| <C N and 
W(X; r) depends on X polynomially, whereas p is ex- 
ponential in X. Therefore we expanded p(X + r) « 
p(X) exp(— rd x s) and replaced u>(x — r/7V; r) — » w(x; r). 

Equation (|1]) reduces the problem of the quasi- 
stationary probability distribution to the problem of clas- 
sical dynamics of an auxiliary Hamiltonian system with 
equations of motion 

x = J2 r rw ( x; r ) ePr ' p = J2 r 5xu; ( x; r ) ( ePr - x ) -( 5 ) 



The distribution p(X) is determined by the mechani- 
cal action of the auxiliary system s(x). In the quasi- 
stationary regime s = 0, i.e., H = in Eq. We notice 
that J2 T w ( x ; r ) r d x s < for i?(x, <9 x s) = 0, and therefore 
s(x) decreases if the point x shifts along a mean-field tra- 
jectory © [provided d x s ^ 0] [l9l. [20|. Since the mean- 
field trajectories go to xa, the action s(x) is minimal at 
x^. Respectively, p(X) is maximal for X^ = iVx^, as 
expected on physical grounds. 

In the spirit of the method of optimal fluctuation 
13, the extinction rate exponent is determined by 
the minimum of s(x) on the extinction hyperplane. From 
the above arguments, the minimum is reached at the ex- 
tinction state X£. Therefore 

/oo 
dtp±. (6) 
-oo 

Equation ([B]) corresponds to the intuitive picture in which 
the most probable fluctuation leading to extinction starts 
from the stable state and brings the system to the extinc- 
tion state, cf. Fig. [TJ The respective optimal Hamilto- 
nian trajectory, Eq. (O, goes from the Hamiltonian fixed 
point (x^,p = 0) to the fixed point (xs,p,g). 

We now find the final momentum p,s. Since action s(x) 
is maximal with respect to x^e at xg, (ps)i^E = 0. To 
find (jps)e we note that, if w(x;r) smoothly vary with 
x, then quite generally, from Eq. ([2]) u>(x;r) oc xe for 
te 7^ 0, and from H = 

[x £ 1 «;(x;r)] x ^ xs {cxp[(p s ) £ r £ ]-f} = 0. (7) 

Equation ([7]) has a trivial solution (ps)e = 0. How- 
ever, there are no Hamiltonian trajectories that would go 
from (xa,P = 0) to (xs,p = 0). Indeed, using Eq. ([3]) 
one can show that trajectories that go to (xs,p = 0) lie 
on the manifold xe = Q,Pi=£E = 0. This manifold does 
not contain the point (x^,p = 0). The trajectory that 
gives the exponent Q goes to (xs,ps) with (ps)e ^ 0, 
as found earlier for specific models 0, 0, [HJ . Therefore 
near x,$ the quasi-stationary distribution p as a function 
of x steeply varies with xe, p oc exp[— N(j>s)eXe]- 

The above analysis can be extended to systems with ef- 
fective constraints, which are implicit in functions iu(x; r) 
and, for example, give extra integrals of motion. In this 
case the above conclusions about change; in partic- 
ular it is no longer necessary to have (ps)i=>tE = 0. To 
gain intuition into this change and its dramatic effect on 
Q we will consider the problem of disease extinction in 
the SIS model. In this model the numbers of susceptible 
and infected individuals X\ and X% change because of 
birth and death, with rates 

W(X;(1,0)) =Nfi, W{X;(-l,0))= f tX 1 , 
W{X;(0,-l))=pX 2 , (8) 

and because of infection on contact and recovery, with 
those recovered immediately becoming susceptible 0]. 
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The corresponding rates are 



where X! = (C/i+x)/(,5+/i) with f x = exp[{j3C - x){t- 
ta)]; constants C, t 2 should be found by matching the 



W(X; (-1, 1)) - (3XtX 2 /N, W(X; (1, -1)) - xX 2 .(9) so i ut i ons given by Eqs. {TTJ, d; ■ If we set C = R, 
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Disease extinction occurs where X 2 = AT^ = 0. For 
the infection reproductive rate Rq > 1, where i?o = 
f3/{p + x), the system has an endemic equilibrium = 
X^/iV = (Rn , 1 — Rq 1 ) ■ It coexists with the disease- 
free stationary state xs = ~K$/N = (1,0). 

Much work has been done on the SIS model in the 
limit p = where the total population does not fluctu- 
ate, x\ + X2 = 1 @, B 0, @, E3- Here, the Hamiltonian 
system ([5]) has effectively one degree of freedom. A direct 
substitution shows that on the optimal trajectory p 2 = 
and pi = h\(f3x\ / k) , which gives 



Q^=o = N(\nR -l+R. Q 



(10) 



In this case (ps)e = (ps)2 = whereas (ps)i = lni?o- 
This is in contradiction with the general result for ex- 
tinction in unconstrained systems and is a consequence 
of the conservation of the total population. 

We now consider the situation where the total popula- 
tion is fluctuating, albeit slowly, that is the characteristic 
birth-death rate pi <C x. Still we assume that p ^> W, 
so that the distribution is quasi-stationary. The Hamil- 
tonian trajectory for extinction consists of three almost 
straight sections IT, T2, T3 shown in Fig.[2jb). Sections 
T1,T3 correspond to slow motion characterized by time 
fi , whereas motion in section T2 is fast, with typical 
time (/3 — x) -1 . A direct substitution shows that motion 
in section Tl is described by equations 



Pi = P2 



= ln l + e^-^ 



X2 



,,l>-2 



Rn 



111) 
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then Xi — > R , 2:2 ™* 0, and p\ — > for t — > 00, and 
the trajectory approaches section T3. On section T3 

P2 = -\n(RoXi), xi = l-exp[-fi(t-t 3 )], (13) 

while [pi|,iC2 — ► for \i — > 0. 

The solutions match if at the end of section Tl and 
at the beginning of section T2 we have pi = P2 = 
— (lni?o)/2. At the end of section T3 we have x — > 
xs = (1, 0) and p — > ps = (0, — lni? ), as expected from 
the general analysis of unconstrained extinction problem. 
The extinction rate exponent is 



N{Rl' 2 - 1) 2 /R 



(14) 



This value, which is obtained in the limit (i — ► 0, is 
smaller than Q for fi = 0, cf. Eq. (JTOJ) . The disconti- 
nuity with respect to n shows the fragility of the result 
obtained by disregarding fluctuations of the total popu- 
lation. 

In Fig.[3]we compare the values of Q obtained for p, = 
and for fj, — > 0. Also shown are the results of numerical 
simulations obtained for (i = and for small nonzero 
p. They are in excellent agreement with the analyti- 
cal results. As illustrated in Fig. [2ja), for Wt <C 1 the 
probability distribution accumulates linearly in time near 
x\ = X\/N = 1 in the extinction plane, X2 = 0. Away 
from the extinction plane, for discrete X2 > 1, the dis- 
tribution is quasi-stationary. For small X2 it has a peak 

— 1/2 

along x\ at « i? where the asymptotic extinction 
path approaches the plane X2 = 0. 



while \xi — R 1 | < \i \t\ in Eq. (fTTj) is arbitrary]. 
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FIG. 2: (Color online), (a) A snapshot of the probability 
p(X) near the extinction plane X2 = for the SIS model; p is 
quasi-continuous in X\/N. The data of simulations refer to 
pt = 9, 7?o = 4, p = p/(p + /t) = 0.1. For t = the system 
was at Xyi, the total number of particles was iV = 50. (b) 
Asymptotic optimal Hamiltonian trajectories for extinction 
for p, — > (solid line) and p — (dashed line). 

Motion in section T2 can be described by setting p = 
in Eqs. This gives 

p 2 = InC, pi = ln(Ci? xi), xi + x 2 = C, (12) 
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FIG. 3: (Color online). The switching exponent Q for the 
SIS model of epidemics. The solid and dashed lines show the 
results for p — > [Eq. {T4}] and p = [Eq. (JTDJ], respectively. 
The data points are obtained from the numerical solution of 
the master equation for the total initial populations N = 50 
and ./V = 100, which made it possible to directly extract the 
exponent Q. 

The fragility in the SIS model results from the inap- 
plicability of a perturbation theory in the fluctuations of 
the total population. We now show that a perturbation 
theory quite generally breaks down in the problem of ex- 
tinction in systems with effective constraints (integrals of 
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motion); the fragility in the SIS model follows from this 
analysis. We assume that, because of the constraint, on 
the optimal extinction trajectory (ps)i^E ^ at least 
for one i. A perturbation changes the elementary transi- 
tion rates inEq. JTJ, py(X;r) W(K; r) + fj,W^ (X; r), 
with /i <C 1 for a small perturbation. The Hamiltonian 
in the eikonal approximation for the probability distribu- 
tion, Eq. (fj|, is respectively modified, H — > H + /iifW. 
To first order in /z the resulting change of the extinc- 
tion exponent can be calculated along the trajectory 
x(t),p(t) of the unperturbed Hamiltonian [2l| . 

/oo 
dtH^(K(t),p{t)) (15) 
-OO 

if (1) (x, p) = V (x; r) [cxp(pr) - 1] . 

Because p(t) exponentially decays for t — > — oo (where 
x — > x^), the integral over time in Eq. (|15p does not 
diverge at the lower limit . 

In order for the perturbation not to destroy the ex- 
tinction state, must satisfy condition j2]), and then 
quite generally w^ 1 ' = W^'/N oc xe for xe — > and 
te i= 0. Since on the optimal extinction trajectory XE^t) 
exponentially decays for t — > oo, the integral (fT5j) of the 
terms with te ^ in fi^ 1 ) converges on the upper limit. 
However, because by assumption at the endpoint of the 
optimal trajectory Pi^e for some i, and because gen- 
erally 

w (1) (x 5 ;r) ^ for r E = 0, (16) 



the Hamiltonian if^ 1 ) remains nonzero for t — > oo and 
overall the integral Eq. (fT5")) diverges. 

The divergence of Q' 1 ' means that the optimal extinc- 
tion trajectory changes nonperturbatively, as does also 
the rate exponent Q, i.e., the extinction rate is fragile 
with respect to the corresponding perturbation. Popu- 
lation fluctuations in the SIS model provide an example 
of such a perturbation, as seen from the comparison of 
Eqs. |(5J) and (fTS)) . We note that the divergence of the 
perturbation theory does not emerge in the problem of 
switching over a saddle point, since p — > as the optimal 
trajectory approaches the saddle point pol |. 

In conclusion, we have considered the exponent in the 
extinction rate Q and demonstrated that it may be frag- 
ile. A small perturbation (oc fj,) can change it signifi- 
cantly, Q for /x — ► differs from Q for /j, = 0. The fragility 
is related to the discontinuity of the quasi-stationary ex- 
tinction current and the related steep slope of the qua- 
sistationary probability distribution near the extinction 
state. A formal condition for the onset of fragility is de- 
rived. Explicit results are obtained for the broadly used 
SIS model of epidemics, and it is shown that this model 
is fragile with respect to fluctuations of the total popula- 
tion. The analytical results are quantitatively confirmed 
by simulations. 
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